# marten et al - 2019 - run experiments.R
#
# runs the simulations for the paper "Exploring the General Equilibrium Costs of
# Sector-Specific Environmental Regulations" by Alex Marten, Richard Garbaccio,
# and Ann Wolverton
#
# this script should be run from the top level directory
#



# ------------------------------------------------------------------------------
# set options
# ------------------------------------------------------------------------------

# benchmark file
benchmark_file = "default_aggregation.gdx"

# mapping file associated with the benchmark
map.file = "default_aggregation.gms"

# shock size - billions of benchmark year dollars
shock.size = 0.1

# affected sources in sensitivity analyses
sensitivity.vintage = "all"

# shock type for sensitivity analyses
sensitivity.shock = "np.air"

# type of image format in which to save plots
plot.type = "tiff"

# dpi of saved plots
plot.dpi = 800

# major grid line color
grid.line.color = "grey75"

# location of the R script's directory relative to the sage model
proj.directory = file.path(".")

# colors for plots
sage.palette = c("#9C9F84","#A97D5D","#F7DCB4","#5C755E","#6D929B","#B7AFA3",
                 "#FF9966")

# plot size for output files
plot.width = 8
plot.height = 4.5



# ------------------------------------------------------------------------------
# load libraries
# ------------------------------------------------------------------------------

# load in the R utilities for sage
source(file.path("utilities","R_utilities.R"))

library(ggplot2)



# ------------------------------------------------------------------------------
# set global variables
# ------------------------------------------------------------------------------

# get the list of sectors and years in sage
sectors = get.sage.sectors(map.file)



# ------------------------------------------------------------------------------
# check for necessary support directories
# ------------------------------------------------------------------------------

# create the images directory if it doesn't exist
if (!dir.exists(file.path(proj.directory,"images")))
  dir.create(file.path(proj.directory,"images"))


# create the results directory if it doesn't exist
if (!dir.exists(file.path(proj.directory,"results")))
  dir.create(file.path(proj.directory,"results"))



# ------------------------------------------------------------------------------
# function definitions
# ------------------------------------------------------------------------------


write.productivity.shock = function(shock.type,shock.size,sector,ex_share=-1) {

  # determine the appropriate flag for the shock type
  if (shock.type=="hicks")
    type = 1
  else if (shock.type=="np.air")
    type = 2
  else if (shock.type=="capital")
    type = 3
  else if (shock.type=="labor")
    type = 4

  # create a gams input file containing the policy shock
  fc = file(file.path("examples","shock_definition.gms"),open="wt")
  writeLines(paste0('s_reg(\"',sector,'\") = YES;'),con=fc)
  writeLines(paste0('ex_shock_share = ',ex_share,';'),con=fc)
  writeLines(paste0('shock_type = ',type,';'),con=fc)
  writeLines(paste0('shock_size = ',shock.size,';'),con=fc)
  writeLines('',con=fc)
  close(fc)

}



# this function loops through each sector writes a policy file for the requested
# shock, runs the model, and stores the welfare change, which is then returned
run.experiments = function(vintage,shock.type,options="",shock.size=0.1,
                           benchmark_file="default_aggregation.gdx") {

  # determine share of shock that applies to extant capital based production
  # from the vintage flag
  if (vintage=="all")
    ex_share = -1
  else if (vintage=="new")
    ex_share = 0
  else if (vintage=="extant")
    ex_share = 1

  # name of output file with the model results
  output_file = "cost_experiments.csv"

  # print a header denoting the start of a round of experiments
  print("")
  print(paste0("running ",shock.type," for ",vintage," capital..."))

  # initialize container for the results
  results = NULL

  # loop over each sector in the model and shock it
  for (i in 1:length(sectors)) {

    # current sector
    sector = sectors[i]

    # write the policy file for the sector
    sector.shock.size = write.productivity.shock(shock.type,
                                                 shock.size,
                                                 sector,
                                                 ex_share=ex_share)

    # run the policy scenario
    # note that "sector" is not a real option for the model, if there is an
    # error in the solve r prints out a warning message with the system call
    # and this way we can identify which sector had the error
    print(paste0("running ge model for ",sector,"..."))
    run.model(benchmark_file=benchmark_file,
              output_file=output_file,
              policy_file="examples/productivity_shock.gms",
              parameter_file=parameter_file,
              gdx_baseline_file=gdx_baseline_file,
              gdx_save=0,
              options=paste(options," --sector=",sector))

    # pause to allow gams to fully release the csv file with the results
    Sys.sleep(1.5)

    # load the results
    ge.policy = read.csv(file.path("output",output_file),stringsAsFactors=FALSE)

    # extract the welfare change
    temp = ge.policy[ge.policy$parameter=="ev",c("region","household","value")]
    temp$sector = sector

    results = rbind(results,temp)

  }

  # add the shock name to the results
  results$shock = shock

  # add a readable list of capital type effected to the results
  results$vintage = vintage

  return(results)

}



# divide the ge welfare costs by the ex ante cost - this is not straight forward
# since it needs to take into account changes in extant vs. new production and
# regional differences in production and interest rates
get.ge.cost.ratio = function(results,bau,shock.size=0.1) {

  # aggregate the results across regions and households to get the national
  # welfare change for each sector's run
  ratio = aggregate(results$value,
                    by=list(results$sector,results$shock,results$vintage),
                    sum)
  names(ratio) = c("sector","shock","vintage","value")


  for (i in 1:nrow(ratio)) {

    # get baseline output levels by captial type, region, and year
    y_ex = bau$value[bau$sector==ratio$sector[i] & bau$parameter=="y_ex"]
    y = bau$value[bau$sector==ratio$sector[i] & bau$parameter=="y"]

    # get baseline returns on new capital by region and year
    pr = bau$value[bau$parameter=="pr"]

    # determine share of shock that applies to extant capital based production
    # from the vintage flag
    if (ratio$vintage[i]=="all")
      ex_share = NA
    else if (ratio$vintage[i]=="new")
      ex_share = 0
    else if (ratio$vintage[i]=="extant")
      ex_share = 1

    # get the vector of years in the data -> the length of this is number of
    # years times the number of regions
    times = bau$time[bau$parameter=="pr"]

    # if the shock affected all vintages set the extant share to the region and
    # year specific share of sectoral output from extant capital production
    if (any(is.na(ex_share)))
      ex_share = sum(y_ex[times==min(times)])/
                 sum(y_ex[times==min(times)]+y[times==min(times)])

    # for each year calculated the present value of the ex ante shock and add
    # it to the running total for the sector
    if (ratio$vintage[i]=="all")
      sector.shock.size = sum(shock.size/sum(y[times==min(times)]+
                                             y_ex[times==min(times)])*y_ex*pr,
                            na.rm=TRUE)+
                        sum(shock.size/sum(y[times==min(times)]+
                                           y_ex[times==min(times)])*y*pr,
                            na.rm=TRUE)
    else
      sector.shock.size = sum(ex_share*shock.size/sum(y[times==min(times)]+
                                                      y_ex[times==min(times)])*
                                                  y_ex*pr,
                              na.rm=TRUE)+
                          sum((1-ex_share)*shock.size/sum(y[times==min(times)]+
                                                          y_ex[times==min(times)])*
                                                      y*pr,
                              na.rm=TRUE)

    # save the ratio of the ge costs to the ex ante shock size
    ratio$value[i] = ratio$value[i]/sector.shock.size

  }

  return(ratio)

}



# ------------------------------------------------------------------------------
# run the ge experiments
# ------------------------------------------------------------------------------

# run the dynamic baseline
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv")
bau.dynamic = read.csv(file.path("output","baseline_dynamic.csv"),
                       stringsAsFactors=FALSE)

# sage parameter file
parameter_file = "parameters.gms"

# sage baseline results file
gdx_baseline_file = "baseline_dynamic.gdx"

# container for the results
results = NULL

# combinations of effected capital types to run
capital.types = c("all","new","extant")

# shock types to run
shock.types = list("hicks","capital","labor","np.air")

# run the experiments
for (shock in shock.types) {
  for (capital in capital.types) {
    results = rbind(results, run.experiments(capital,
                                             shock,shock.size=shock.size,
                                             benchmark_file=benchmark_file))
  }
}



# ------------------------------------------------------------------------------
# function for plotting results
# ------------------------------------------------------------------------------

# this is the ordering of sectoral results in the plot
man.sectors = c("bom","cem","chm","con","cpu","fbm","fmm","pmm","prm","tem","wpm")
eng.sectors = c("col","cru","ele","gas","ref")
oth.sectors = c("agf","min","trn","ttn","wsu")
sector.order = rev(c(man.sectors,eng.sectors,oth.sectors))

# function to plot the results
plot.results = function(data,x.limits=NA) {

  # only keep the sectors to be plotted
  data = data[data$sector %in% sector.order,]
  data$sector = factor(data$sector,sector.order)

  # resource extraction sectors don't have extant capital so extant only shocks
  # have zero cost, therefore we remove those experiments from the plot
  data$value[data$value==-100] = NA

  # creat data for the min and max
  bars = data.frame(sector=unique(data$sector))
  for (i in 1:nrow(bars)) {
    bars$ymin[i] = min(data$value[data$sector==bars$sector[i]],na.rm=TRUE)
    bars$ymax[i] = max(data$value[data$sector==bars$sector[i]],na.rm=TRUE)
  }

  if(!any(is.na(x.limits))) {
    y.max = x.limits[2]-4
  } else {
    y.max = max(bars$ymax)
  }

  # plot the results
  p = ggplot()+
    geom_point(data=data,aes(x=sector,y=value,color=scenario,shape=scenario),
               stat="identity",size=3,stroke=2)+
    geom_hline(aes(yintercept=0))+
    geom_line(aes(y=rep(round(y.max)+2,length(man.sectors)),
                  x=length(sector.order)-(1:length(man.sectors))+1))+
    geom_line(aes(y=rep(round(y.max)+2,length(eng.sectors)),
                  x=length(sector.order)-length(man.sectors)-(1:length(eng.sectors))+1))+
    geom_line(aes(y=rep(round(max(y.max))+2,length(oth.sectors)),
                  x=1:length(oth.sectors)))+
    geom_text(aes(y=round(y.max)+4,
                  x=sector.order[length(sector.order)-round(length(man.sectors)/2)],
                  label="Manufacturing",angle=270))+
    geom_text(aes(y=round(y.max)+4,
                  x=sector.order[length(sector.order)-length(man.sectors)-round(length(eng.sectors)/2)],
                  label="Energy",angle=270))+
    geom_text(aes(y=round(y.max)+4,
                  x=sector.order[round(length(oth.sectors)/2)+1],
                  label="Other",angle=270))+
    coord_flip()+
    scale_color_manual(values=sage.palette)+
    scale_shape_manual(values=c(0,1,2,5))+
    labs(x="Sector",y="Percent Difference Between GE Costs and Engineering Costs")+
    guides(color=guide_legend(title="Scenario",override.aes=list(shape=c(0,1,2,5))),
           shape="none")+
    theme(axis.line.x      = element_line(colour="black"),
          axis.line.y      = element_line(colour="black"),
          panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_line(color=grid.line.color),
          panel.grid.minor = element_blank(),
          panel.border     = element_blank(),
          panel.background = element_blank(),
          legend.box.just  = "left",
          legend.key       = element_blank(),
          legend.background = element_rect(fill="transparent",color="transparent"))

  # add optional limits for the x axis
  if (!any(is.na(x.limits)))
    p = p+scale_y_continuous(limits=x.limits)

  print(p)

}


# ------------------------------------------------------------------------------
# plot the input composition sensitivities
# ------------------------------------------------------------------------------

data = get.ge.cost.ratio(results,bau.dynamic,shock.size=shock.size)
data$value = (data$value-1)*100

data = data[data$vintage==sensitivity.vintage,]

# rename the input composition scenarios
names(data)[names(data)=="shock"] = "scenario"
data$scenario[data$scenario=="hicks"] = "Neutral"
data$scenario[data$scenario=="capital"] = "Capital Only"
data$scenario[data$scenario=="labor"] = "Labor Only"
data$scenario[data$scenario=="np.air"] = "Nestor & Pasurka"

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Nestor & Pasurka",
                                       "Neutral",
                                       "Capital Only",
                                       "Labor Only"))

plot.results(data,x.limits=c(0,37))

ggsave(file.path("images",paste0("figure_6a.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=plot.width,height=plot.height)



# ------------------------------------------------------------------------------
# plot the vintage sensitivities
# ------------------------------------------------------------------------------

data = get.ge.cost.ratio(results,bau.dynamic,shock.size=shock.size)
data$value = (data$value-1)*100

data = data[data$shock==sensitivity.shock,]

# rename the vintage scenarios
names(data)[names(data)=="vintage"] = "scenario"
data$scenario[data$scenario=="all"] = "All Sources"
data$scenario[data$scenario=="extant"] = "Existing Sources"
data$scenario[data$scenario=="new"] = "New Sources"

# order the experiments for the plots
data$scenario = factor(data$scenario,c("All Sources",
                                       "Existing Sources",
                                       "New Sources"))

data$value[data$sector %in% c("gas","col","agf","cru","min")] = NA

plot.results(data,x.limits=c(0,37))

ggsave(file.path("images",paste0("vintage_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# create table with mean ge effects
# ------------------------------------------------------------------------------

# get the list of shocks and vintages runs
shocks = unique(results$shock)
vintages = unique(results$vintage)

# initialize a space for the table
table = matrix("",nrow=length(shocks)+1,ncol=length(vintages)+1)
table[1,1] = "Input Shares"
for (i in 2:nrow(table))
  table[i,1] = shocks[i-1]
for (j in 2:ncol(table))
  table[1,j] = vintages[j-1]

# poopulate the table with means and std. devs. in parantheses
for (i in 2:nrow(table)) {
  for (j in 2:ncol(table)) {
    values = results[results$shock==shocks[i-1] & results$vintage==vintages[j-1] & results$sector %in% sector.order,]
    values = get.ge.cost.ratio(values,bau.dynamic,shock.size=shock.size)
    values$value = (values$value-1)*100
    values$value[is.infinite(values$value)] = NA
    table[i,j] = paste0(round(mean(values$value,na.rm=TRUE)),"% (",
                        round(sqrt(var(values$value,na.rm=TRUE))),")")
  }
}

# save the table
write.csv(table,file.path(proj.directory,"results","ge_costs.csv"))



# ------------------------------------------------------------------------------
# run the shock size sensitivities
# ------------------------------------------------------------------------------

# values for the shocks to be included in the sensitivity analysis
shock.sizes = c(.4,.8,1.2,1.6,2)

# container for the results intialized with the default results
results.size = results[results$shock==sensitivity.shock &
                         results$vintage==sensitivity.vintage,]
results.size = get.ge.cost.ratio(results.size,bau.dynamic,shock.size=shock.size)
results.size$scenario = shock.size

for (s.size in shock.sizes) {

  # run the experiments
  temp = run.experiments(sensitivity.vintage,sensitivity.shock,shock.size=s.size)

  # convert to percentage difference from engineering costs
  temp = get.ge.cost.ratio(temp,bau.dynamic,shock.size=s.size)
  temp$scenario = s.size

  # save the results
  results.size = rbind(results.size,temp)

}



# ------------------------------------------------------------------------------
# plot the shock size sensitivity results
# ------------------------------------------------------------------------------

data = results.size

# convert to an index relative to the ratio at the default shock size
for (s in sectors)
  data$value[data$sector==s] = data$value[data$sector==s]/
                               data$value[data$sector==s & data$scenario==shock.size]

# convert shock size from billion $ to million $
data$scenario = data$scenario*1000

temp = aggregate(data$value[data$sector %in% man.sectors],
                 by=list(data$scenario[data$sector %in% man.sectors]),mean)
names(temp) = c("scenario","value")
temp$Sector = "Manufacturing"
plot.data = temp

temp = aggregate(data$value[data$sector %in% eng.sectors],
                 by=list(data$scenario[data$sector %in% eng.sectors]),mean)
names(temp) = c("scenario","value")
temp$Sector = "Energy"
plot.data = rbind(plot.data,temp)

temp = aggregate(data$value[data$sector %in% oth.sectors],
                 by=list(data$scenario[data$sector %in% oth.sectors]),mean)
names(temp) = c("scenario","value")
temp$Sector = "Other"
plot.data = rbind(plot.data,temp)

plot.data$Sector = factor(plot.data$Sector,c("Manufacturing","Energy","Other"))

ggplot(plot.data)+
  geom_line(aes(x=scenario,y=value,color=Sector),size=1.5)+
  scale_color_brewer(palette="Pastel2")+
  labs(x="Annual Ex Ante Engineering Costs [million 2016$]",
       y=paste0("Percent Difference Between GE Costs and Engineering\n Costs Relative to $",
                shock.size*1000," Million Policy"))+
  scale_y_continuous(limits=c(.9,1))+
  theme(axis.line.x      = element_line(colour="black"),
        axis.line.y      = element_line(colour="black"),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color=grid.line.color),
        panel.grid.minor = element_blank(),
        panel.border     = element_blank(),
        panel.background = element_blank(),
        legend.box.just  = "left",
        legend.key       = element_blank(),
        legend.background = element_rect(fill="transparent",color="transparent"))

ggsave(file.path("images",paste0("size_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# run the labor supply sensitivities
# ------------------------------------------------------------------------------

# container for the results
results.l = results[results$shock==sensitivity.shock &
                      results$vintage==sensitivity.vintage,]
results.l = get.ge.cost.ratio(results.l,bau.dynamic,shock.size=shock.size)
results.l$scenario = "Default"

# sage baseline results file
gdx_baseline_file = "baseline_dynamic.gdx"

# create a sage parameter file without labor-leisure choice
# done by copying the default sage parameters file and zeroing out (essentially)
# the calibrated labor supply elasticities
parameter_file = "parameters_no_ll.gms"
file.copy(file.path("data","parameters.gms"),
          file.path("data","parameters_no_ll.gms"),overwrite=TRUE)
fc = file(file.path("data","parameters_no_ll.gms"),"at")
writeLines("lse_sub = 0.000001;",con=fc)
writeLines("lse_inc = -0.000001;",con=fc)
close(fc)

# run the dynamic baseline without labor-leisure choice
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file)
bau.dynamic.l = read.csv(file.path("output","baseline_dynamic.csv"),
                         stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.l,shock.size=shock.size)
temp$scenario = "Perfectly Inelastic Labor Supply"

# save the results
results.l = rbind(results.l,temp)

# create a sage parameter file with a larger labor supply elasticity
# done by copying the default sage parameters file and adding larger calibrated
# labor supply elasticities
parameter_file = "parameters_hi_lse.gms"
file.copy(file.path("data","parameters.gms"),
          file.path("data","parameters_hi_lse.gms"),overwrite=TRUE)
fc = file(file.path("data","parameters_hi_lse.gms"),"at")
writeLines("lse_sub = 0.4;",con=fc)
close(fc)

# run the dynamic baseline with more elastic labor supply
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file)
bau.dynamic.l = read.csv(file.path("output","baseline_dynamic.csv"),
                         stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.l,shock.size=shock.size)
temp$scenario = "High Labor Supply Elasticity"

# save the results
results.l = rbind(results.l,temp)



# ------------------------------------------------------------------------------
# plot the labor supply sensitivities
# ------------------------------------------------------------------------------

data = results.l
data$value = (data$value-1)*100

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Default",
                                       "Perfectly Inelastic Labor Supply",
                                       "High Labor Supply Elasticity"))

plot.results(data)

ggsave(file.path("images",paste0("l_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# run the value-added substitution elasticity sensitivities
# ------------------------------------------------------------------------------

# container for the results
results.kl = results[results$shock==sensitivity.shock &
                       results$vintage==sensitivity.vintage,]
results.kl = get.ge.cost.ratio(results.kl,bau.dynamic,shock.size=shock.size)
results.kl$scenario = "Default"

# sage baseline results file
gdx_baseline_file = "baseline_dynamic.gdx"

# create a sage parameter file with lower se_kl
parameter_file = "parameters_se_kl.gms"
file.copy(file.path("data","parameters.gms"),
          file.path("data","parameters_se_kl.gms"),overwrite=TRUE)
fc = file(file.path("data","parameters_se_kl.gms"),"at")
writeLines("se_kl(s) = .5*se_kl(s);",con=fc)
close(fc)

# run the dynamic baseline with lower se_kl
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file)
bau.dynamic.kl = read.csv(file.path("output","baseline_dynamic.csv"),
                          stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.kl,shock.size=shock.size)
temp$scenario = "Low KL Elasticity"

# save the results
results.kl = rbind(results.kl,temp)

# create a sage parameter file with higher se_kl
file.copy(file.path("data","parameters.gms"),
          file.path("data","parameters_se_kl.gms"),overwrite=TRUE)
fc = file(file.path("data","parameters_se_kl.gms"),"at")
writeLines("se_kl(s) = 1.5*se_kl(s);",con=fc)
close(fc)

# run the dynamic baseline with lower se_kl
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file)
bau.dynamic.kl = read.csv(file.path("output","baseline_dynamic.csv"),
                          stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.kl,shock.size=shock.size)
temp$scenario = "High KL Elasticity"

# save the results
results.kl = rbind(results.kl,temp)

# create a sage parameter file with higher se_kl
file.copy(file.path("data","parameters.gms"),
          file.path("data","parameters_se_kl.gms"),overwrite=TRUE)
fc = file(file.path("data","parameters_se_kl.gms"),"at")
writeLines("se_kl(s) = 1.0001;",con=fc)
close(fc)

# run the dynamic baseline with lower se_kl
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file)
bau.dynamic.kl = read.csv(file.path("output","baseline_dynamic.csv"),
                          stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.kl,shock.size=shock.size)
temp$scenario = "Unit KL Elasticity"

# save the results
results.kl = rbind(results.kl,temp)



# ------------------------------------------------------------------------------
# plot the se_kl sensitivities
# ------------------------------------------------------------------------------

data = results.kl
data$value = (data$value-1)*100

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Default",
                                       "Low KL Elasticity",
                                       "High KL Elasticity",
                                       "Unit KL Elasticity"))

plot.results(data)

ggsave(file.path("images",paste0("se_kl_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# run the k structure sensitivities
# ------------------------------------------------------------------------------

# container for the results
results.k = results[results$shock==sensitivity.shock &
                      results$vintage==sensitivity.vintage,]
results.k = get.ge.cost.ratio(results.k,bau.dynamic,shock.size=shock.size)
results.k$scenario = "Default"

# sage baseline results file
gdx_baseline_file = "baseline_dynamic.gdx"

# default sage parameter file
parameter_file = "parameters.gms"

# run the dynamic baseline without putty-clay capital
run.model(benchmark_file=benchmark_file,output_file="baseline_dynamic.csv",
          parameter_file=parameter_file,options=" --putty_clay=0")
bau.dynamic.k = read.csv(file.path("output","baseline_dynamic.csv"),
                          stringsAsFactors=FALSE)

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size,options=" --putty_clay=0")

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.dynamic.k,shock.size=shock.size)
temp$scenario = "No Putty-Clay Capital"

# save the results
results.k = rbind(results.k,temp)



# ------------------------------------------------------------------------------
# plot the k structure sensitivities
# ------------------------------------------------------------------------------

data = results.k
data$value = (data$value-1)*100

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Default",
                                       "No Putty-Clay Capital"))

plot.results(data)

ggsave(file.path("images",paste0("k_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# run the tax sensitivities
# ------------------------------------------------------------------------------

benchmark_files = c("default_aggregation_no_taxes.gdx",
                    "default_aggregation_no_tl.gdx",
                    "default_aggregation_no_tk.gdx")

# container for the results intialized with the default results
results.tax = results[results$shock==sensitivity.shock &
                        results$vintage==sensitivity.vintage,]
results.tax = get.ge.cost.ratio(results.tax,bau.dynamic,shock.size=shock.size)
results.tax$scenario = benchmark_file

# sage baseline results file
gdx_baseline_file = "baseline_dynamic.gdx"

# default sage parameter file
parameter_file = "parameters.gms"

for (b_file in benchmark_files) {

  # run the dynamic baseline
  run.model(benchmark_file=b_file,output_file="baseline_dynamic.csv")
  bau.dynamic.tax = read.csv(file.path("output","baseline_dynamic.csv"),
                             stringsAsFactors=FALSE)

  # run the experiments
  temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                         benchmark_file=b_file,shock.size=shock.size)

  # convert to percentage difference from engineering costs
  temp = get.ge.cost.ratio(temp,bau.dynamic.tax,shock.size=shock.size)
  temp$scenario = b_file

  # save the results
  results.tax = rbind(results.tax,temp)

}



# ------------------------------------------------------------------------------
# plot the tax sensitivity
# ------------------------------------------------------------------------------

data = results.tax
data$value = (data$value-1)*100

data$scenario[data$scenario=="default_aggregation.gdx"] = "Default"
data$scenario[data$scenario=="default_aggregation_no_taxes.gdx"] = "No Taxes"
data$scenario[data$scenario=="default_aggregation_no_tl.gdx"] = "No Labor Tax"
data$scenario[data$scenario=="default_aggregation_no_tk.gdx"] = "No Capital Tax"

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Default",
                                       "No Taxes",
                                       "No Labor Tax",
                                       "No Capital Tax"))

plot.results(data)

ggsave(file.path("images",paste0("tax_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)



# ------------------------------------------------------------------------------
# run static vs. dynamic comparison
# ------------------------------------------------------------------------------

# container for the results intialized with the default results
results.static = results[results$shock==sensitivity.shock &
                           results$vintage==sensitivity.vintage,]
results.static = get.ge.cost.ratio(results.static,bau.dynamic,shock.size=shock.size)
results.static$scenario = "Dynamic"

# sage parameter file
parameter_file = "parameters_static.gms"

# run the static baseline
run.model(benchmark_file=benchmark_file,output_file="baseline_static.csv",
          parameter_file=parameter_file)
bau.static = read.csv(file.path("output","baseline_static.csv"),
                             stringsAsFactors=FALSE)

# sage baseline results file
gdx_baseline_file = "baseline_static.gdx"

# run the experiments
temp = run.experiments(sensitivity.vintage,sensitivity.shock,
                       shock.size=shock.size)

# convert to percentage difference from engineering costs
temp = get.ge.cost.ratio(temp,bau.static,shock.size=shock.size)
temp$scenario = "Static"

# save the results
results.static = rbind(results.static,temp)



# ------------------------------------------------------------------------------
# plot the static vs. dynamic comparison
# ------------------------------------------------------------------------------

data = results.static
data$value = (data$value-1)*100

# order the experiments for the plots
data$scenario = factor(data$scenario,c("Dynamic","Static"))

plot.results(data)

ggsave(file.path("images",paste0("static_sensitivity.",plot.type)),
       device=plot.type,dpi=plot.dpi,width=8,height=4.5)
